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Abstract 



Reversibility, weak reversibility and deficiency, detailed and complex balancing are generally 
not "encoded" in the kinetic differential equations but they are realization properties that may 
imply local or even global asymptotic stability of the underlying reaction kinetic system when 
further conditions are also fulfilled. In this paper, efficient numerical procedures are given for 
finding complex balanced or detailed balanced realizations of mass action type chemical reaction 
networks or kinetic dynamical systems in the framework of linear programming. The procedures 
are illustrated on numerical examples. 

Keywords: reaction kinetic systems, mass action kinetics, linear programming 
AMS classification: 80A30 chemical kinetics 

1 Introduction 

Chemical reaction networks form a wide and intensively studied class of nonnegative polynomial 
systems [39]. Beside describing pure chemical reactions, CRNs are often used to model complex 
biological mechanisms [30], or models of application fields seemingly far from chemistry such as 
mechanical or electrical systems [31]. The increasing and extended interest for this field is shown 
by the fact that numerous surveys and tutorials have been published even in journals where the 
primary scope is not chemistry [34, 4, 1]. 

A special subclass of CRNs is the set of reaction networks where the so-called mass action law 
(MAL) is assumed to be valid. The strongest theoretical results concerning the dynamical properties 
of CRNs are available for MAL kinetics. Chemical Reaction Network Theory (CRNT) initiated by 
the results published e.g. in the pioneering papers and lectures [24, 15, 16, 14] deals with the 
relations between structures/parameters and dynamical behaviour of CRNs. Ever since, there has 
been continuous development in this field with significant contributions such as in [6, 7, 19, 33]. 

The principle of detailed balance has been applied in thermodynamics for long, which says that 
the forward and backward rates of a certain molecular process must be equal at the equilibrium [3]. 
For chemical reactions, this means that the forward and reverse rates of each reversible elementary 
reaction step should be equal at the equilibrium state [28]. In the 1970s, a special set of CRNs called 
complex balanced systems were identified that exhibit "regular" behaviour of the concentration 
trajectories in the sense that at least local stability of the positive equilibrium points can be proved 
and certain "exotic" dynamical phenomena such as periodic or chaotic solutions can be ruled out 
[23, 13]. 

The widely known fundamental dogma of chemical kinetics says that parametrically and/or 
structurally different CRNs can give rise to the same mathematical model in ODE form [32]. In 
other words, seemingly very different CRNs may have exactly the same function in terms of the 
time evolution of the species concentrations. Such networks will be called dynamically equivalent 
in this paper. In [35] the notions of dense and sparse realizations were introduced for dynamically 
equivalent reaction networks containing the maximal or minimal number of reactions with a fixed 
set of complexes. Moreover, a computational method in the framework of mixed integer linear 
programming (MILP) was proposed to compute such realizations. In [36], important properties of 
dense realizations were described: it was shown that the structure (i.e. the unweighted directed 
graph) of a dense realization of a given CRN is unique, and the unweighted directed reaction 
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graph of any other realization is the subgraph of the dense realization if the set of complexes is 
fixed. Furthermore, additional MILP based methods were given for computing reversible CRN 
realizations, and realizations with the minimal number of complexes. 

Since, similarly to the number of linkage classes, reversibility, weak reversibility and deficiency, 
detailed and complex balancing are generally also not "encoded" in the kinetic differential equa- 
tions but they are realization properties, it is of interest whether a kinetic system or a CRN with 
multiple possible structures has a detailed/complex balanced realization or not. The existence of 
such realization(s) may imply local or even global asymptotic stability of the underlying reaction 
kinetic system when further conditions are also fulfilled. 

The outline of the paper is the following. In section 2, the basic notions in the fields of basic 
CRN properties, kinetic realizability of polynomial systems, and linear programming (LP) will be 
summarized. Section 3 contains the main contributions where the computation of CRN realizations 
with detailed balance and complex balance is casted into LP problems. Section 4 contains illustrative 
examples, and finally, the most important conclusions can be found in section 5. 

2 Preliminaries 

2.1 Nonnegative systems 

The concepts and basic results in this subsection are mostly taken from [4]. A function / = 
[fi ■ ■ ■ fn] T '• [0, oo) n —¥ R n is called essentially nonnegative if, for all i = 1, . . . , n, fi(x) > for 
all x G [0, oo) n , whenever X{ = 0. In the linear case, when f(x) = Ax, the necessary and sufficient 
condition for essential nonnegativity is that the off-diagonal entries of A are nonnegative (such a 
matrix is often called a Metzler-matrix) . 

Consider an autonomous nonlinear system 

x = f(x), x(0) = x (1) 

where / : X — > R n is locally Lipschitz, X is an open subset of R n and xq € X . Suppose that 
the nonnegative orthant [0, oo) n = R™ C X. Then the nonnegative orthant is invariant for the 
dynamics (1) if and only if / is essentially nonnegative. 

2.2 Mass action reaction networks 

The overview in this subsection is largely based on [35] and [36]. A CRN obeying the mass action 
law is a closed system where chemical species Xj, i = 1, n take part in r chemical reactions. The 
concentrations of the species denoted by Xi, (i = 1, n) form the state vector, i.e. Xj = [Xj]. The 
elementary reaction steps have the following form: 

n n 

->-^/%Xj, j = l,...,r (2) 

i=l i=l 

where ajj is the so-called stoichiometric coefficient of component X, in the jth reaction, and fin is 
the stoichiometric coefficient of the product X^. The linear combinations of the species in eq. (2), 
namely a ij^-i an d Y^t=i Pij^i f° r J = 1) • • • > r are called the complexes and are denoted by 
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Ci, C2, . . . , C m . Note that i/ie stoichiometric coefficients are always nonnegative integers in classical 
reaction kinetic systems. The reaction rates of the individual reactions can be described as 

n n 

I'M) /••,n^ X -" 'oil'-'/ > i = 1 >-,r (3) 

i=l i=l 

where fcj > is the reaction rate coefficient of the jth reaction. The reaction rate and reaction rate 
coefficient of the reaction Cj — > Cj will also be denoted by Pij(x) and kij, respectively, whenever 
it is more convenient. 

In our computations, the following form will be used for the description of the dynamics of 
CRNs obeying the mass action law [14]: 

± = Y-A k - ^(x) (4) 

where x = [X] G W 1 is the concentration vector of the species, Y € R nxm stores the stoichiometric 
composition of the complexes, A k 6 W nxm contains information about the structure of the reaction 
network, and ip : R n 1— > R m is a monomial- type vector mapping given by 

n 

^j( x ) = II x ' QlJ ' i = 1 '---' m ( 5 ) 

i=l 

where ctij = \Y]ij- To denote the right hand side of eq. (4), the notation M = Y ■ A k will also be 
used. The structure of Y and A k is the following. The ith column of Y contains the composition 
of complex Cj, i.e. Y^ is the stoichiometric coefficient of Cj corresponding to the specie Xj. Aj~ 
is a column conservation matrix (i.e. the sum of the elements in each column is zero), called the 
Kirchhoff matrix of the CRN, defined as 

It is important to note that the pair (Y,Ak) uniquely characterizes a particular realization of a 
given CRN. The invariance of the nonnegative orthant for the dynamics of CRNs is shown by 
several authors (see, e.g. [39, 4]). 

To handle the exchange of materials between the environment and the reaction network, the so- 
called "zero-complex" can be introduced and used which is a special complex where all stoichiometric 
coefficients are zero i.e., it is represented by a zero vector in the Y matrix [14]. 

For the reaction Cj — > Cj, the corresponding reaction vector is defined as 

v k = [Y]. tj - [Yy (7) 

where [Y]. t denotes the ith column of Y. Similarly to reaction rate coefficients, whenever it is more 
practical, Vij denotes the reaction vector corresponding to the reaction Cj — > Cj. 

The rank of a reaction network denoted by s is defined as the rank of the vector set H = 
{v\,V2 ■ ■ ■ ,v r } where r is the number of reactions. The elements of H span the so-called stoichio- 
metric subspace denoted by S, i.e. S = spanj^i, V2 ■ ■ ■ , v r }. The positive stoichiometric compatibility 
class containing a concentration xq is the following set [15]: 

(x + 5) n Rl 
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where denotes the positive orthant in M. n . 

The deficiency d of a reaction network is defined as [14, 15] 

d = m — I — s (8) 

where m is the number of complexes in the network, I is the number of linkage classes (graph 
components) and s is the rank of the reaction network. 

Similarly to [14] and many other authors, the following weighted directed graph (often called 
Feinberg- Horn- Jacks on graph) is assigned to the reaction network (2). The directed graph D = 
(Vd, Ed) of a reaction network consists of a finite nonempty set Vd of vertices and a finite set of 
ordered pairs of distinct vertices called directed edges. The vertices correspond to the complexes, i.e. 
Vd = {Ci, C2, . . . C m }, while the directed edges represent the reactions, i.e. (Cj, Cj) £ Ed if complex 
Cj is transformed to Cj in the reaction network. The reaction rate coefficients kj for j = 1, . . . ,r 
in (3) are assigned as positive weights to the corresponding directed edges in the graph. A reaction 
network is called reversible, if whenever it contains the reaction Cj — > Cj , then the reverse reaction 
Cj < — Cj is also present in the CRN. A reaction network is called weakly reversible, if each complex 
in the reaction graph lies on at least one directed cycle (i.e. if complex Cj is reachable from complex 
Cj on a directed path in the reaction graph, then Cj is reachable from Cj on a directed path). An 
important point of the well-known Deficiency Zero Theorem [15] says that the ODEs of a weakly 
reversible deficiency zero CRN are globally stable with a known logarithmic Lyapunov function for 
all positive values of the reaction rate coefficients. 

If the corresponding forward and backward reaction pairs form a cycle in the directed graph of 
a reversible CRN, then this cycle of reversible edge-pairs will be called a circuit (to distinguish it 
from a directed cycle). E.g. the CRN in Fig. 2 c) contains a circuit of length 3, while the CRNs in 
Figs. 2 a) and b) do not contain any circuit. 

2.3 Kinetic realizability of nonnegative polynomial systems 

An autonomous polynomial nonlinear system of the form (1) is called kinetically realizable or simply 
kinetic, if a mass action reaction mechanism given by eq. (4) can be associated to it that exactly 
realizes its dynamics, i.e. f(x) = Y ■ A). ■ ip{x) where matrix Y has nonnegative integer elements, 
tp contains the monomials according to (5), and A}~ is a valid Kirchhoff matrix. In such a case, 
the pair {Y,A^) will be called a realization of the kinetic system (1). As it is expected from linear 
algebra, the same polynomial system may have many parametrically and/or structurally different 
realizations, and this is particularly true for kinetic systems coming from application domains other 
than chemistry. The problem of kinetic realizability of polynomial vector fields was first examined 
and solved in [25] where the constructive proof contains a realization algorithm that produces the 
weighted directed graph of a possible associated mass action mechanism. According to [25], the 
necessary and sufficient condition for kinetic realizability is that all coordinates functions of the 
right hand side of (1) must have the form 

fi(x) = -Xigi(x) + hi(x), i = l,...,n (9) 

where gi and hi are polynomials with nonnegative coefficients. Roughly speaking, condition (9) 
means that kinetic systems cannot contain negative cross-effects. From this, it is easy to see that all 
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nonnegative linear systems are kinetic, since a linear system characterized by a Metzler matrix where 
only the diagonal elements can be negative is obviously in the form of (9). Moreover, classical Lotka- 
Volterra (LV) systems that are known to be nonnegative, are always kinetic, too. However, there 
are many essentially nonnegative polynomial systems that are not directly kinetic, since some of the 
monomials in /j that do not contain X{ have negative coefficients. To circumvent this problem, one 
possible way is to embed the differential equations into the LV class that can be done algorithmically 
[22]. However, this solution usually results in the significant dimension increase of the state space 
(i.e. in the increase of species in the corresponding kinetic realization). If the equations and the 
possible initial conditions guarantee that all state variables remain strictly positive throughout 
the solutions, then a more advantageous method is a simple state-dependent time-rescaling of the 
equations that preserves the important qualitative properties of the system and the solutions [37, 21]. 
The above summary clearly shows that the class of systems that are kinetic or can be transformed 
to kinetic form is really wide within nonnegative systems. 

The very short description of the realization algorithm presented in [25] is the following. Let us 
write the polynomial coordinate functions of the right hand side of a kinetic system (1) as 

Ti n 

f i (x)=J2™i j l[x b3k , i = l,...,n (10) 

3=1 k=l 

where 7*j is the number of monomial terms in /j. Let us denote the transpose of the iih standard 
basis vector in M. n as and let Bj = [bji . . . bj n ]. Then, the necessary steps (transcribed into an 
easily implementable form) for constructing a realization are given by: 

Algorithm 1 from [25] 

For each i = 1, . . . , n and for each j = 1, . . . , do: 

1. Cj = Bj + sign(mij) • 

2. Add the following reaction to the graph of the realization 

n n 
k=l k=l 

with reaction rate coefficient \rriij\, where Cj = [cji ... Cj n ]. 



As we can see later in section 4, Algorithm 1 often produces a non-minimal reaction structure with 
more reactions and/or complexes than the minimal number needed for the kinetic representation of 
the studied polynomial system. 

2.4 Complex balance and detailed balance 

Definition 2.1. A CRN realization (Y, Aj.) is called complex balanced at x* € if 

A k i;(x*) = 0. (11) 







Definition 2.2. A reversible CRN realization (Y,Ak) is called detailed balanced at x* £ M™, if 

Pij(x*) = pji(x*), such that C{ ^ Cj exists (12) 

i. e. the forward and reverse reaction rates for each reversible reaction are equal at x* . 

Clearly, if a CRN is complex balanced at x*, then x* is an equilibrium point, i.e. YA) s ip(x*) = 0. 
Furthermore, it is easy to check that detailed balancing at x* implies complex balancing at x* for 
reversible networks, but the converse is generally not true. 

Definition 2.3. A reversible CRN realization (Y,Ak) is called complex balanced (detailed balanced), 
if it is complex balanced (detailed balanced) at each positive steady state. 

It is important to remark, that both the complex balance and detailed balance properties depend 
on the structure of the reaction digraph and on the numerical values of the reaction rate coefficients 
[5]. In [17], necessary and sufficient conditions were given for detailed balancing (see also [29] for a 
clear explanation of conditions) using the circuit and spanning forest conditions below. The circuit 
conditions require that the product of the reaction rate coefficients should be equal taken in either 
directions along a circuit. The spanning forest conditions say that 

H(%) 7 - =i[(k ji y« (13) 

where the products are taken for each reaction step in the spanning forest of the reaction digraph 
and 7 denotes the solutions of the equation 

lijVjj = 0; ^ or an hJ sucn that Ci — > Cj exists, (14) 

and Vij is the reaction vector of the reaction Cj — > Cj. 

It is worth summarizing the following important properties and relations of detailed and complex 
balancing collected from [13, 23, 17, 20, 5, 10]. 

PI If there is no circuit in the reaction graph, then the spanning forest condition is alone necessary 
and sufficient for detailed balancing. 

P2 If the deficiency of the CRN is zero, then the circuit condition is necessary and sufficient for 
detailed balancing. 

P3 If the deficiency of the CRN is zero, and there is no circuit in the reaction digraph, then the 
reaction is detailed balanced. 

P4 If the circuit conditions are satisfied in a reversible CRN containing at least one circuit, then 
detailed balancing and complex balancing are equivalent. 

P5 If a CRN is complex balanced (detailed balanced) at any positive x* then it is complex balanced 
(detailed balanced) at all other positive equilibrium points. 

P6 If a CRN is complex balanced then it is weakly reversible. 
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P7 A CRN is complex balanced for any positive values of reaction rate coefficients if and only if a) 
it is weakly reversible, and b) the deficiency of the system is zero. 

P8 If a CRN is complex balanced then there is precisely one equilibrium point in each stoichiometric 
compatibility class. Furthermore, each equilibrium point is at least locally stable with a known 
strict Lyapunov function. 

Additionally, we have to mention the so-called Global Attractor Conjecture for which no counterex- 
ample has been found, but the construction of a complete proof seems to be technically challenging 
[5]. The conjecture says that for any complex balanced CRN and any initial condition x(0), the 
equilibrium point x* is a global attractor in the corresponding positive stoichiometric compatibility 
class. The conjecture further supports the need for a method to find complex balanced realizations 
of kinetic systems, if stability properties of such systems are to be studied. The circuit and spanning 
tree conditions for detailed balancing are hard to insert into an optimization framework in their 
original form, but another simple algebraic condition will allow the use of linear programming for 
this purpose. Moreover, property P5 will enable us to use an arbitrary positive steady state to find 
a detailed balanced or complex balanced realization if it exists. 

2.5 Linear programming 

Linear programming (LP) is an optimization technique, where a linear objective function is mini- 
mized/maximized subject to linear equality and inequality constraints [8, 9]. LP has been widely 
used in chemistry and chemical engineering in the fields of system analysis [2, 26], simulation [18] 
and design [27]. The standard form of linear programming problems that will be used in this paper 
is the following 



where y € M™ is the vector of decision variables, c S M n A S W xn , b € M p are known vectors and 
matrices, and '>' in (17) means elementwise nonstrict inequality. 

The feasibility of the simple LP problem (15)-(17) can be checked using the following necessary 
and sufficient condition [8]. 

Theorem 2.1. Consider the auxiliary LP problem 



minimize c y 



(15) 



subject to: 
Ay = b 



(16) 
(17) 



v 




(18) 



8=1 



subject to: 
Ay + z = b 

y>0,z>0 



(19) 
(20) 
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where z GMP is a vector of auxiliary variables. Let us assume (without the loss of generality) that 
b > 0. There exists a feasible solution for the LP problem (15) -(17) if and only if the auxiliary LP 
problem (18) -(20) has optimal value with Z{ = for i = 1, . . . ,p. 

It is easy to see that y = 0, z = b is always a feasible solution for (18)-(20). The above theorem 
will be useful for establishing that no detailed or complex balanced realization exists in certain 
cases. 



3 Finding balanced realizations using linear programming 



In this section, we assume that an initial CRN realization (Y^\A^) or a kinetic polynomial system 
is given together with an arbitrary positive steady state x* that has been determined analytically 
or simply through simulations. If we start from a kinetic polynomial system, we use Algorithm 
1 for generating an initial realization (Y^\ Aj}'). In any case, 

M (i) = y(i) . A W Let us denote 

tp{x*) simply by ip*. 

3.1 Constraints representing mass action dynamics 

Similarly to [35], let us introduce the following notations. The Kirchhoff-matrix of the searched 
realization is written as 



,(*) 



-a lx 

021 



Ol2 
-0t22 



aim 

0<2m 



a-ml «m2 • • • 

Then, the optimization variable, y is defined columnwise as 
[ an a 2 i a 3 i . 



(21) 



y 



«ml 0,12 



(m—l)m a mm 



(22) 



Taking into consideration the diagonal elements in (21) with a negative sign enables us to constrain 
all elements of y to be nonnegative. Let [W] i . and [W], ■ denote the zth row and jth column of an 
arbitrary matrix W, respectively. Let us define the following matrices (with known entries) as 



Y l 



[Y«]. fl 
1 



[y(i 



,i-l 



[yd)] [yd)] 



[yd)] 
1 



i(n+l)xra 



(23) 



M 



Md) 

... 



G 



a(n+l)xm 



Then the equality constraints encoding mass action dynamics can be written as 

A e i ■ y = B e i 



(24) 



(25) 
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where 



A 



el 



Y 1 ... 
Y 2 ... 



im(n+l)xm 2 



(26) 



... Y m 
(with zeros representing zero blocks of size (n + 1) x m), and 



(27) 



[Ml 

The inequalities expressing sign constraints for the reaction rate coefficients are simply 

y>o. 



(28) 



The choice of the objective function In principle, the objective function (15) can be any 
linear function of y that is determined by the choice of the vector c. In section 4, we will minimize 
or maximize the sum of reaction rate coefficients. The minimization of the sum corresponds to the 
choice of c = ^[1 ... 1] T in the standard LP-problem (15)-(17), while the maximization is obtained 
by selecting c = — 1[1 ... 1] T . Note that the minimum of the sum of reaction rate coefficients is 
bounded if the LP constraints are feasible, but the maximum is not necessarily, and this must be 
kept in mind. 

3.2 Additional constraints for complex balancing 

Using the definition of a complex balanced steady state, the equality constraints for complex bal- 
ancing are 

-an^ + a 12 ^* 2 + ■■■ + a im V4 = (29) 
a 2 lV'* - «22'02 H 1" a2mV4 = ( 30 ) 



amllpl H H aWtm-l^m-l ~ ^mmV'm = 

Clearly, these constraints can be written as 

A e2 ■ y = B e2 , 

where 



.4 



c2 



ipl 



m—l 



im-1 



and denotes a 1 x k row vector of zeros. 



m—2 



rpl O 171 - 1 ipZ ... 



m—l 



B e2 = G 



(31) 
(32) 

(33) 



10 



3.3 Further constraints for detailed balancing 

The most suitable form of detailed balancing constraints is taken from [38]. According to this, a 
given steady state x* is detailed balancing if and only if 



G ■ AT = A k ■ G 



(34) 



where G = diag(?/>*). It is easy to see that (34), if satisfied, implies reversibility of the obtained 
reaction network. Eq. (34) encodes a maximum of m ^ ^ independent equations that can be 
written into the linear programming problem as 



since 



,4 



1p*y{i-l)m+j - ^jy(j-l)m+i = 0) V * > 3, 
V(j-i)m+i according to eq. (22). 



(35) 



4 Examples 

4.1 Multiple detailed balanced realizations of a simple irreversible network 




Figure 1: Simple deficiency 1 reaction network 

Consider the simple reaction network shown in Fig. 1 with parameters k\ = 1, hi = 1.5. This 
CRN was also used as an example in [36] in a different context. The matrices characterizing complex 
composition and graph structure are 



Y 



3 2 
3 1 



Thus, 

M = Y ■ A k 
The corresponding ODE model is 





" -1 







1 


-1.5 







1.5 


3 


-1.5 


" 


-3 


1.5 






(36) 



X\ — 1.53/1 



X2 



+ 1.5s? 



(37) 



(38) 
(39) 
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It's easy to compute that e.g. x* = [1.6725 1.3275] r is a steady state for the system (38)-(39). The 
corresponding values of the monomials are tp* = ip(x*) = [2.3393 4.6786 3.7134] T . Clearly, the 
realization shown in Fig. 1 cannot be complex balanced or detailed balanced, since it is not weakly 
reversible. 

In the first step, let us define the objective function as 



1 9 



(40) 



i=l 



which is the Ll-norm of the reaction rate coefficient vector (note that the diagonal elements of 
with a negative sign contain the sum of the reaction rate coefficients in the corresponding column). 
Minimizing the objective function h in (40) subject to the conditions (25), (28) and (35) gives the 
detailed balance solution 



Vopt 



= [ 1 1 0.5 


0.5 








Lg realization 










" -1 


0.5 





Y' = Y, A' k = 


1 


-0.5 

















The detailed balance condition at the steady state can be checked as 

1 . ( x *f = 0.5 • (x^) 3 = 2.3393. 



(41) 



(42) 



(43) 



which is a reversible, deficiency realization. Therefore, by applying the Deficiency zero theorem, 
the dynamics of the initial irreversible network is also globally stable with a known entropy-like 
Lyapunov function. 

In the second step, let us maximize h in (40). The obtained realization is now given by 



Y" = Y, A'L 



-1.5 0.9449 
-1.5 1.8899 
1.5 1.5 -2.8348 



(44) 



The deficiency of the network in this case is 1. The detailed balance property is fulfilled in this case 
as 



1.5 • {x* 2 f = 0.9449 • {x\fxl = 3.5088, 
1.5 • {x\f = 1.8899 • {x\) 2 x* 2 = 7.0179. 



(45) 
(46) 



The reaction digraphs of the above computed detailed balanced realizations can be seen in Fig. 2 
a) and b), respectively. 

The motivation for trying to minimize and maximize the LI norm of the reaction rate coefficients 
to get potentially different detailed balanced realizations came from [12] and [11]. We remark that 
the exact conditions under which the sparsest solution of an underdetermined set of linear equations 
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3X„ 



0.5 



a) 




1.8899 



2X 1+ X 



b) 



0.9333 

► OV 

0.4667 1 




2X 1+ X 2 



C) 



Figure 2: Dynamically equivalent detailed balanced realizations of the CRN shown in Fig. 1. a) 
deficiency zero realization; b), c) deficiency one realizations 



is the minimal LI norm solution are not fulfilled in our case, but the two obtained solutions are 
structurally different, and our purpose here was only to illustrate that different detailed balanced 
realizations of the same (very simple) kinetic system may exist. 

It's worth mentioning that the detailed balancing constraints described in section 3.3 can easily 
be combined with the computation of dense and sparse realizations (see [35, 36]). However, the 
problem in this case becomes NP-hard. A dense detailed balanced realization can be seen in Fig. 2 
c) that is determined by the following matrices: 



Y 



y, k 



-1.0333 0.4667 0.0630 
0.9333 -0.5667 0.1260 
0.1000 0.1000 -0.1890 



The detailed balancing condition is met in this case, too, since 



0.063 
0.126 



x\) 2 x* 2 = 0.2339 



0.9333 • (x* 2 y = 0.4667 • (x 
0.1 • (x*f 
0.1 • (xlf 



*\3 



2.1833 



0.4679 



(47) 



(48) 
(49) 
(50) 



The average running time of the LP based methods (using the linprog command) was 0.06 second, 
while the MILP based algorithm (using the freely available GLPK solver) found the dense realization 
in 1.04 second in the MATLAB computation environment on a notebook computer with an 1.66 
GHz Intel Atom N280 Processor. 



4.2 Finding complex balanced realization of a kinetic polynomial system 

The following polynomial system is given: 

X\ = x\ — X\X2 + X3X4 — 2X1X2X3 
X2 = x\ — X\X2 + 2X3X4 — \X\x\x3 

±3 = —2x\ + x\X2 — X\x\x3 + 2x 4 (51) 
£4 = X1X2 — X3X4 + Ax\x\x3 — 3^4 

It can be seen that (51) is essentially nonnegative and kinetic. After running Algorithm 1, we 
obtain an initial network with 19 complexes and 16 reactions that is visible in Fig. 3, where the 
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numbering of complexes (in rectangles) is the following: 



1 : 2X 3 , 2 : X 3 + X 4 , 3 : Xi + 2X 3 , 4 : X 2 + 2X 3 , 

5 : X 3 , 6 : Xi + X 3 + X 4 , 7 : X 2 + X 3 + X 4 , 

8 : Xi + X 2 , 9 : Xi + 2X 2 + X 3 , 10 : X 1; 11 : X 2 , 

12 : Xi + X 2 + X 4 , 13 : Xi + X 2 + X 3 , 14 : 2X 2 + X 3 , 

15 : X x + 2X 2 , 16 : Xx + 2X 2 + X 3 + X 4 , 

17 : 3X 4 , 18 : X 3 + 3X 4 , 19 : 2X 4 

The objective function that was minimized was also the sum of reaction rate coefficients, i.e. 



(52) 
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Figure 3: Reaction network realizing eq. (51) obtained using Algorithm 1. Only those reaction 
rate coefficients are indicated that are different from 1. 



361 

%) = 5 I> ( 53 ) 
i=l 

The constraints were given by eqs. (25), (28) and (32). The joint running time of Algorithm 1 and 
the solution of the LP problem with 361 variables was 0.45 second on the same hardware/software 
environment as in the previous example. The obtained complex balanced and thus weakly reversible 
realization of the initial CRN is visible in Fig. 4. The isolated complexes (i.e. the ones with no 
incoming and no outgoing directed edges) are naturally omitted from the realization. It can be 
checked that the deficiency of the complex balanced realization is (in sharp contrast with the 
deficiency of 12 of the initial network in Fig. 3). Therefore, the autonomous system (51) has well- 
characterizable equilibrium points in the positive orthant that are globally stable with a known 
entropy-like Lyapunov function [15]. 

An attempt to find a detailed balanced realization using the constraints described in section 
3.3 is unsuccessful in this case. After checking the feasibility condition given in Theorem 2.1, we 
can deduce that such detailed balanced realization does not exist with the complex set listed in eq. 
(52). 
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3X — 2X,^ X 1+ X ; 

\ / 
X +2X +X^X +X 

1 2 3 3 4 



Figure 4: Complex balanced realization of the CRN shown in Fig. 3. All reaction rate coefficients 
are 1. 

5 Conclusions 

Polynomial-time methods were proposed in this paper for the computation of complex balanced 
and detailed balanced realizations of kinetic systems. The input of the methods is a given kinetic 
polynomial system (or an initial CRN) together with an arbitrary steady state point. The algorithms 
are based on standard linear programming where the characteristics of mass action dynamics are 
represented in the form of linear constraints. If solvable, the methods provide sufficient conditions 
for the existence of reversible or weakly reversible realizations in the cases of detailed balance or 
complex balance, respectively. The existence of such realizations usually guarantee strong stability 
properties for the underlying ODEs [5]. If a kinetic polynomial system is given, then a key step of 
the proposed methods is the algorithm described in [25] for generating an initial realization. The 
methods can be easily combined with the computation of dense and sparse realizations [36, 35]. 
The efficiency of the methods was demonstrated on numerical examples. 
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